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ABSTRACT 

PDF  methods  provide  a  means  of  calculating  the  properties  of  turbulent  reacting  flows.  They  have  been 
successfully  applied  to  many  turbulent  flames,  including  some  with  finite-rate  kinetic  effects.  Here  the 
methods  are  reviewed  with  an  emphasis  on  computational  issues  and  on  their  application  to  turbulent 
combustion. 


1 .  INTRODUCTION 

PDF  methods  have  been  applied  to  a  variety  of  turbulent  flows  both  with,  and  without,  combustion. 
Generally,  single-phase,  low-Mach-number  flows  have  been  considered,  in  which  radiation  is  not  a  major 
factor.  For  such  flows,  the  fundamental  dependent  variables  are  the  velocities  U(x,t)  and  the  compositions 

ct»(x,t)  (e.g.,  the  species  mass  fractions  and  enthalpy).  Different  probabilistic  approaches  to  modelling 

turbulent  flows  can  be  categorized  according  to  the  statistics  of  U(x,t)  and  $(x,t)  that  are  considered.  For 

example,  in  a  mean-flow  closure  <H(x,t)>  and  <<Ji(x,t)>  are  the  primary  dependent  variables:  in  second- 

t  t  i 

order  closures  the  variances  <uiuj>,  <<t)a<t)p>  and  covariances  <Ui<|)a>  are  also  included.  Angled  brackets 

f 

denote  means  (i.e.,  mathematical  expectations)  and  u,(x,t)  =  Ui(x,t)-<Ui(x,t)>  and  <t>a(x,t)  =  <|>a(x,t) 
-«t>a(x.t)>  are  the  fluctuating  components  of  Ui  and  <|>a,  respectively. 

In  pdf  methods,  the  dependent  variable  is  a  probability  density  function  (pdf)  or  joint  pdf  of  U(x,t), 

<i>(x,t).  The  pdf  contains  information  equivalent  to  all  the  moments;  and  hence  pdf  methods  are,  in  this  sense, 
more  comprehensive  than  moment  closures  (e.g.  second-order  closures).  The  methods  that  have  proved 
most  successful  are  based  on  one-point,  one-time  pdfs,  which  contain  information  at  each  point  in  the  flow 
separately,  but  no  joint  information  at  two  or  more  distinct  points. 

In  the  last  15  years,  pdf  methods  have  advanced  from  being  only  of  theoretical  interest  to  a  small  group 
of  specialists,  to  being  a  practical  approach  for  calculating  the  properties  of  turbulent  reactive  flows.  As  the 
review  below  indicates,  in  additional  to  having  been  applied  to  idealized  flames  and  simple  laboratory  flames, 
the  methods  have  been  applied  to  flames  requiring  multistep  chemical  kinetics  (e.g.  Pope  1981a,  Chen  & 
Kollmann  1988a),  as  well  as  to  computationally  difficult  flows  (e.g.  that  in  the  cylinder  of  a  spark-ignition 
engine,  Haworth  &  El  Tahry  1989a, b). 

The  purpose  of  this  paper  is  to  review  the  work  on  pdf  methods,  with  some  emphasis  on  the  numerical 
issues,  and  on  the  applications  to  turbulent  combustion. 

In  the  next  section  the  different  pdf  methods  are  described,  along  with  the  associated  modelling.  Monte 
Carlo  methods  have  proved  to  be  the  most  successful  means  of  solving  pdf  transport  equations.  The  essence 
of  these  solution  techniques  is  described  in  Section  3.  The  theoretical  foundations  of  pdf  methods  (including 
the  modelling  and  Monte  Carlo  solution  algorithms)  are  comprehensively  described  by  Pope  (1985).  The 


purpose  of  Sections  2  and  3  is  to  describe  briefly  the  principal  features,  with  no  attempt  at  rigor.  Section  4 
reviews  the  applications  of  pdf  methods  to  turbulent  diffusion  flames  and  premixed  flames.  (Recent 
applications  to  constant-density  inert  flows  have  been  reviewed  by  Pope  1989a.)  In  the  Discussion  and 
Conclusions  some  of  the  outstanding  problems  and  future  directions  are  assessed. 

2 .  PDF  METHODS 


2.1  Definitions  and  Properties 

Let  <|>  denote  the  value  of  a  composition  variable  (the  mass  fraction  of  oxygen,  for  example)  at  a  particular 
location  Xo  anti  time  to  in  a  turbulent  reactive  flow.  It  is  supposed  that  the  flow  can  be  realized  any  number 
of  times,  and  the  time  t  is  measured  from  the  initiation  of  the  flow.  Thus  from  each  realization  we  obtain  a 

value  of  <|);  and,  given  the  nature  of  turbulence,  these  values  are,  in  all  probability,  different.  In  other  words 

<|>  is  a  random  variable.  On  a  given  realization,  it  is  not  possible  to  determine  beforehand  the  value  of  <(>  that 
will  obtain.  But  it  is  possible  to  ascribe  probabilities  to  its  value  being  in  a  given  interval:  this  can  be  done 
through  the  pdf. 

For  every  random  variable  we  introduce  an  independent  (sample-space)  variable:  in  particular,  \|/  is  the 
sample-space  variable  corresponding  to  (j).  The  cumulative  distribution  function  (cdf),  F^y),  is  then  defined 
as  the  probability  that  <j>  is  less  than  \)/: 


F^y)  =  Prob{<|)<\|/}  •  (1) 

And  the  pdf  of  <|),  f<j)(\) /),  is  defined  by 

f4>(V)sT"  *W)-  (2) 

d\| f 

While  F<j>(\j/)  is  a  probability  function,  f<t>(\|0  is  a  probability  density  function.  That  is,  f<j,(\j/)  is  the 
probability  per  unit  y/  of  the  event  (j)=\|/:  or,  equivalently,  f<j>(\j/)d\j/  is  the  probability  of  the  event 

The  three  fundamental  properties  of  the  pdf  (in  addition  to  Eq.  2)  are: 

f<t»(V)>0,  (3) 

since  probabilities  are  non-negative; 

oo 

J  fyty)  d\j/=l  ,  (4) 

-OO 

since  Prob{<(K<»}  =  1  and  Prob  {<(><-«>}  =  0;  and,  for  any  (non-pathological)  function  Q(<J>): 

OO 

<Q>  =  J  f<|>(¥)  Q(¥)d¥  • 


(5) 


This  last  equation  shows  that  if  the  pdf  is  known,  then  the  expectation  of  any  function  of  the  random  variable 
can  be  calculated.  In  particular  the  mean  «(»  and  the  m-th  central  moment  «j)'m>  (m>l)  can  be  determined 
(if  it  exists). 

For  a  general  turbulent  reactive  flow  we  need  to  consider  a  set  of  a  >  1  composition  variables  <|>  = 
{<t>l»<t>2»  ...,<t>a}-  Accordingly  the  a  sample-space  variables  \j/=  {^1,^2.  —»Vct}  are  introduced,  and  the  / 
joint  pdf  of  &  f$0^),  is  defined  to  be  the  probability  density  of  the  compound  event  <t>  =  ^  (i.e.  $1  =  Vi,  4>2  = 

V2, ...,  <>CT  =  Va)- 

Clearly  the  joint  pdf  defined  at  the  particular  location  Xo  and  time  to,  can  be  defined  at  any  (x,t):  we 
denote  by  f$(l|£;x,t)  the  joint  pdf  of  £(x,t).  It  is  important  to  realize  that  this  is  a  one-point,  one-time  joint 
pdf:  it  contains  no  joint  information  between  $  at  two  or  more  positions  or  times.  The  pdf  method  described 
in  the  next  sub-section  is  based  on  %(\£x,t),  which  is  called  the  composition  joint  pdf. 

The  second  pdf  method  described  (in  Section  2.3)  is  based  on  the  velocity-composition  joint  pdf, 
f(V,\jjr;x,t).  Here,  Y  =  { Vi,V2,V3)  are  the  three  independent  velocity  variables;  and  f  is  the  probability 
density  of  the  compound  event  UI(x,t)  =  Y,  <J)(x,t)  =  \j£}. 

In  the  treatment  of  variable-density  flows,  two  further  probability  functions  prove  useful,  and  are  now 
defined.  By  assumption  (see  Pope  1985),  the  composition  variables  are  sufficient  to  determine  the  fluid 

density.  If  the  composition  is  <fe,  the  density  is  given  by  the  function  pa(4>),  which  can  be  determined  from  a 
thermodynamic  calculation.  Consequently  at  (x,t)  the  fluid  density  is 

p(x,t)  =  pa(<Mx,t]) .  (6) 

Having  made  the  distinction  between  the  different  functions  p(x,t)  and  pCT(&),  we  now  follow  conventional 
(if  imprecise)  notation  and  denote  both  by  p. 

Favre,  or  density- weighted,  pdfs  are  defined  by,  for  example: 

f<l>(V)  =  p(v)f<t>(v)/<p>  •  (?) 

It  then  follows  that  density-weighted  means  are  given  by 

OO 

Q=<pQ>/<p>=  I  Q(\|0?(\i/)d\|/ ,  (8) 


(cf.  Eq.  5).  The  mass  density  function  F  is  defined  by 

F( Y,U£,x;t)  =  p(H£)f(Y,\j£;x,t) .  (9) 


The  use  of  these  functions  is  made  apparent  in  the  next  two  subsections. 


2.2  Composition  Joint  PDF  Equation 


Dopazo  and  O'Brien  (1974)  were  the  first  to  consider  the  transport  equation  for  %(\j/;x,t).  Since  then  a 
number  of  different  derivations  have  been  given  (e.g.  Pope  1976,  Janicka  et  al.  1978a, b,  O’Brien  1980, 
Pope  1985).  Here  we  state  the  result,  and  refer  the  reader  to  Pope  (1985)  for  a  detailed  derivation. 

The  compositions  evolve  according  to  the  conservation  equation 


Dt 


“  T~ ’  +  Sa  , 
P  oxi 


(10) 


where  Ia  is  the  (molecular)  diffusive  flux  of  <t>a,  and  Sa  —  a  known  function  of  $  —  is  the  rate  of  creation 
of  <|)a  due  to  chemical  reaction.  The  pdf  transport  equation  corresponding  to  Eq.  (10)  is 
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8t  9xi 


+  [Sa(U£)4]  = 

8Va 


-  —  7-  [<uj V4]  +  r^“  [<“  T"  '^41  • 

<p>  dxj  1  o\|ta  p  dxi 
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On  the  left-hand  side,  the  first  two  terms  represent  the  rate  of  change  following  the  mean  flow.  The  third 
term  is  —  in  composition  space  —  the  divergence  of  the  flux  of  probability  due  to  reaction.  It  is  the  form  of 

this  term  that  gives  this  pdf  method  the  advantage  over  other  statistical  approaches:  since  £(\jf)  is  known,  7$ 

is  the  subject  of  the  equation,  and  \|/a  is  an  independent  variable,  the  term  contains  no  unknowns.  Thus 
however  complicated  and  non-linear  the  reaction  scheme,  in  the  composition  joint  pdf  equation  the  effect  of 
reaction  is  in  closed  form,  requiring  no  modelling. 

In  contrast  the  terms  on  the  right-hand  side  require  modelling.  The  quantity  <u"l\j£>  is  the  mean  of  the 

Favre  velocity  fluctuation  (u"  =  U-tJ)  conditional  upon  the  event  (J>  =  \(/.  It  represents  the  transport  of  7^  in 
physical  space  by  the  fluctuating  velocity.  Although  there  have  been  other  suggestions,  this  term  is  generally 
modelled  by  gradient  diffusion: 


<p><u'.'  I  \£>7$  =  -  Tj  —  7^_, 

OXj 


(12) 


where  Tt  is  a  turbulent  diffusivity.  Such  gradient  transport  models  are,  of  course,  subject  to  many 
objections,  especially  when  applied  to  variable-density  reactive  flows. 

The  final  term  in  Eq.  (11)  represents  the  effect  of  molecular  mixing.  It  is  generally  treated  by  a  stochastic 
mixing  model  (see  e.g.  Flagan  &  Appleton  1974,  Pope  1982,  1985).  While  some  aspects  of  this  modelling 
are  discussed  below,  the  cited  references  should  be  consulted  for  a  full  account. 

The  composition  joint  pdf  equation  is  not  a  self-contained  model:  mean  momentum  equations  must  be 
solved  for  0;  and  a  turbulence  model  (k-e,  say)  is  needed  to  determine  both  Tt  (~  k2/e)  and  the  mixing  rate 
(~  e/k)  used  in  the  stochastic  mixing  model. 


2.2  Composition  Joint  PDF  Equation 


Dopazo  and  O’Brien  (1974)  were  the  first  to  consider  the  transport  equation  for  %(\j/;x,t).  Since  then  a 
number  of  different  derivations  have  been  given  (e.g.  Pope  1976,  Janicka  et  al.  1978a,b,  O'Brien  1980, 
Pope  1985).  Here  we  state  the  result,  and  refer  the  reader  to  Pope  (1985)  for  a  detailed  derivation. 

The  compositions  evolve  according  to  the  conservation  equation 


Dfoq 
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where  Ia  is  the  (molecular)  diffusive  flux  of  <t>a,  and  Sa  —  a  known  function  of  $  —  is  the  rate  of  creation 
of  <[)a  due  to  chemical  reaction.  The  pdf  transport  equation  corresponding  to  Eq.  (10)  is 
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On  the  left-hand  side,  the  first  two  terms  represent  the  rate  of  change  following  the  mean  flow.  The  third 
term  is  —  in  composition  space  —  the  divergence  of  the  flux  of  probability  due  to  reaction.  It  is  the  form  of 

this  term  that  gives  this  pdf  method  the  advantage  over  other  statistical  approaches:  since  &(\j£)  is  known,  7$ 

is  the  subject  of  the  equation,  and  \|/a  is  an  independent  variable,  the  term  contains  no  unknowns.  Thus 
however  complicated  and  non-linear  the  reaction  scheme,  in  the  composition  joint  pdf  equation  the  effect  of 
reaction  is  in  closed  form,  requiring  no  modelling. 

In  contrast  the  terms  on  the  right-hand  side  require  modelling.  The  quantity  <u"l\|£>  is  the  mean  of  the 

Favre  velocity  fluctuation  (u"  =  U-U)  conditional  upon  the  event  (J)  =  \j/.  It  represents  the  transport  of  7^  in 
physical  space  by  the  fluctuating  velocity.  Although  there  have  been  other  suggestions,  this  term  is  generally 
modelled  by  gradient  diffusion: 


<p><u"  I  \jn>7^  =  -  Pr  — -  7$_, 

1  dxj 


(12) 


where  Tt  is  a  turbulent  diffusivity.  Such  gradient  transport  models  are,  of  course,  subject  to  many 
objections,  especially  when  applied  to  variable-density  reactive  flows. 

The  final  term  in  Eq.  (11)  represents  the  effect  of  molecular  mixing.  It  is  generally  treated  by  a  stochastic 
mixing  model  (see  e.g.  Flagan  &  Appleton  1974,  Pope  1982,  1985).  While  some  aspects  of  this  modelling 
are  discussed  below,  the  cited  references  should  be  consulted  for  a  full  account. 

The  composition  joint  pdf  equation  is  not  a  self-contained  model:  mean  momentum  equations  must  be 
solved  for  0;  and  a  turbulence  model  (k-e,  say)  is  needed  to  determine  both  Tx  (~  k2/e)  and  the  mixing  rate 
(~  e/k)  used  in  the  stochastic  mixing  model. 


=  U+(t)  =U(s+[t],t), 


(15) 


since,  by  definition,  a  fluid  particle  moves  with  the  local  fluid  velocity; 


1 

p(4>+)  9xjJx+ 


(from  Eq.  13);  and. 


(16) 


=  ®a  =  SaCk"*’)  • 


1 

P<ife+) 


(17) 


(fromEq.  10). 

The  connection  between  these  equations  for  the  properties  of  a  fluid  particle,  and  the  equation  for  the 
mass  density  function  F  (Eq.  14)  is  immediately  apparent.  Equation  (14)  can  be  written 


dF 

9t 


+  [F<U+I>]  + 

3xj  ^ 


[F<Ajl>]  +  —  [F<@al>]  =  0, 

3Va 


(18) 


where  the  expectations  are  conditional  on  the  compound  event  (x+(t)  =  x,  U+(t)  =  V,  <|)+(t)  =  \j/}.  Further,  it 
may  be  noticed  that  the  terms  in  braces  in  Eqs.  (16)-(17)  appear  on  the  right-hand  side  of  Eq.  (14)  —  that  is, 
they  need  to  be  modelled  —  whereas  all  other  terms  appear  on  the  left-hand  side  and  are  treated  exactly. 

2.5  Stochastic  Models 

The  standard  approach  to  turbulence  modelling  is  to  construct  constitutive  relations  for  the  unknown 
correlations  (see,  e.g.  Lumley  1978).  In  the  context  of  the  mass  density  function,  this  approach  is  to  model 
the  unknown  conditional  expectations  on  the  right-hand  side  of  Eq.  (14)  in  terms  of  known  quantities,  i.e. 

functions  or  functionals  of  F(V,\j/.,x;t).  But  the  Lagrangian  viewpoint  offers  a  different  approach  to 

modelling:  namely  to  use  stochastic  processes  to  simulate  unknown  contributions  to  U+(t)  and  (j)+(t)  (i.e.  the 
terms  in  braces  in  Eqs.  16  and  17). 

To  illustrate  this  approach  we  consider  LJ*(t)  —  a  stochastic  model  for  U+(t).  If  the  model  is  accurate 
then  U*(t)  is  (statistically)  an  accurate  approximation  to  U+(t).  In  general  the  time  series  U*  is  not 
differentiable.  Consequently  we  express  the  models  in  terms  of  the  infinitesimal  increment 

dil*(t)  =  H*(t+dt)  -  U*(t) ,  (19) 

rather  than  in  terms  of  the  derivative  dU*/dt.  Note  that  for  a  deterministic,  differentiable  process,  (e.g. 
U+(t)),  the  infinitesimal  increment  is  non-random  (i.e.  zero  variance)  and  is  of  order  dt. 

In  view  of  the  equation  for  12+(t)  (Eq.  16),  the  increment  dU*  can  be  written 


dt  +  dU? , 


(20) 


where  (similar  to  U*)  x*  and  <fc*  are  models  of  x+  and  <fc+.  The  stochastic  increment  dUs  models  the  effects 
of  the  fluctuating  pressure  gradient  and  viscous  stresses,  while  the  term  in  dt  is  an  exact  expression  for  the 
effect  of  gravity  and  the  mean  pressure  gradient 

Two  types  of  models  have  been  used.  The  first  type  —  of  which  the  stochastic  mixing  model  is  an 
example  —  are  called  particle  interaction  models.  In  the  terminology  of  stochastic  processes,  these  are  point 
processes.  According  to  these  models,  the  infinitesimal  increment  dHs  is  nearly  always  zero.  But  with 
probability  of  order  dt,  the  increment  is  of  order  unity.  Thus  the  time  series  is  a  piecewise  constant,  with  of 
order  unity  jumps  per  unit  time. 

The  second  type  of  model  use  diffusion  processes  in  which  dUs  is  a  random  variable  with  (conditional) 
mean  and  variance  both  of  order  dt.  Note  that  this  implies  that  the  rms  is  of  order  dt1/2,  and  hence  the 
process  —  though  continuous  —  is  not  differentiable.  The  different  variants  of  the  Langevin  model  are 
diffusion  processes  (see  e.g.  Pope  1983, 1985,  Haworth  &  Pope  1986, 1987a). 

For  more  information  on  this  general  modelling  approach  the  reader  is  referred  to  Pope  (1985),  while  the 
current  status  of  the  Langevin  model  is  described  by  Haworth  &  Pope  (1987a). 

3 .  NUMERICAL  SOLUTION  ALGORITHMS 

The  velocity-composition  joint  pdf  f(V.y:x.t)  is  a  single  function  defined  in  a  multi-dimensional  space. 
In  general  f  depends  on  the  three  velocity  variables,  a  composition  variables,  three  spatial  variables,  and  time 

—  (7+a)  independent  variables  in  all.  In  many  cases  the  dimensionality  may  be  less,  but  still  large.  For 
example,  in  a  statistically  stationary  and  two-dimensional  flow  with  a  single  composition  variable, 

f(Y,yi;xi,X2)  depends  on  6  independent  variables.  The  composition  joint  pdf  f<k(\j£;x,t)  in  general  depends 
on  (4+a)  variables;  but,  for  the  simpler  flow  cited  above,  f<t>(v|ri;xi,x2)  is  a  function  of  just  3  variables. 

Given  the  large  dimensionality  of  joint  pdfs  it  is  clear  that  conventional  grid-based  numerical  methods 
(e.g.  finite  differences)  are  impractical  for  all  but  the  simplest  of  cases.  Just  to  provide  an  accurate 
representation  of  a  function  of  6  independent  variables  is  a  major  task.  Consequently,  although  one  or  two 

finite-difference  solutions  have  been  obtained  for  f<j>(yi;xi,X2)  (e.g.  Janicka  &  Kollmann  1978a,b),  currently 
all  investigators  are  using  Monte  Carlo  methods  instead. 

In  the  next  subsection  the  general  Monte  Carlo  method  devised  by  Pope  (1985)  to  solve  for  the  velocity- 
composition  joint  pdf  is  outlined.  Then,  in  section  3.2,  Monte  Carlo  solution  algorithms  for  the  composition 
joint  pdf  are  reviewed. 

3.1  Monte  Carlo  Method  for  the  Velocity-Composition  Joint  PDF 

The  Monte  Carlo  method  to  solve  the  modelled  equation  for  the  velocity-composition  joint  pdf  is 
conceptually  simple  and  natural.  Rather  than  discretizing  the  space,  we  discretize  the  mass  of  fluid  into  a 
large  number  N  of  representative  or  stochastic  particles.  At  a  given  time  t,  let  M  be  the  total  mass  of  fluid 

within  the  solution  domain.  Then  each  stochastic  particle  represents  a  mass  Am  =  M/N  of  fluid.  The  n-th 

particle  has  position  2^(0,  velocity  LJ(n)(t),  and  composition  <j)(n)(t). 

Starting  from  appropriate  initial  conditions,  the  particle  properties  are  advanced  in  time  by  the  increments 

dx00(t)  =  U(n)(t)  dt ,  (21) 

dU(n>(t)  =  [g-p(#n>)-1V<p>]dt  +  dUs , 


(22) 


djfc(n>(t)  =  S(^(n>)dt  +  d£s , 


(23) 


where  dLJs  and  d<ks  are  the  stochastic  increments  that  simulate  molecular  processes  and  the  fluctuating 
pressure  gradient.  At  symmetry  boundaries  particles  are  reflected;  at  inflow  boundaries  particles  are  added 
with  appropriate  properties;  and,  at  outflow  boundaries  particles  are  discarded.  While  wall  boundaries  have 
been  treated  (Anand  et  al.  1989a)  a  comprehensive  account  of  this  treatment  is  not  available  in  the  literature. 

The  correspondence  between  the  ensemble  of  stochastic  particles  and  the  joint  pdf  has  been  established 
by  Pope  (1985).  The  main  results  are: 

N 

i)  The  expected  density  of  the  stochastic  particles  in  physical  space  (Am  £  <5(x-2c(n)(t)>)  is  equal  to  the 

n=l 

fluid  density  <p(x,t)>. 


ii)  The  joint  pdf  of  the  stochastic  particle  properties  LJ(n)(t),  <t>(n)(t)  is  the  density-weighted  joint  pdf 
7(V,\|fix(n)[t],t). 

iii)  From  particle  properties,  expectations  (e.g.  H(x»t))  can  be  approximated  as  ensemble  averages,  with  a 
statistical  error  of  order  N*1/2. 

Several  implementations  of  this  algorithm,  and  variants  of  it,  have  been  reported.  For  example,  the 
turbulent  jet  diffusion  flame  calculations  reported  in  Section  4  are  performed  using  a  "boundary-layer" 
variant  (see  Pope  1985).  Haworth  &  Pope  (1987b)  report  a  variant  of  the  algorithm  designed  specifically  for 
self-similar  shear  flows.  From  a  numerical  standpoint  this  work  is  of  particular  interest,  because  the 
convergence  of  the  method  (as  N'1/2)  is  demonstrated.  The  basic  algorithm  has  been  implemented  and 
demonstrated  for  statistically  two-dimensional  recirculating  flows  by  Anand  et  al.  (1989b).  Haworth  &  El 
Tahry  (1989a,b)  reports  calculations  of  the  three-dimensional  time-dependent  flow  in  the  cylinder  of  a  spark- 
ignition  engine.  In  this  case  the  pdf  algorithm  is  coupled  to  a  conventional  finite- volume  algorithm  that  is 
used  to  calculate  the  mean  pressure  and  turbulent  time  scale  fields. 

3.2  Monte  Carlo  Algorithms  for  the  Composition  Joint  PDF 

Two  different  algorithms  have  been  proposed  to  solve  the  modelled  transport  equation  for  the 
composition  joint  pdf. 

The  second  of  these  (chronologically)  was  proposed  by  Pope  (1985),  and  is  similar  to  that  described 
above  for  the  velocity-composition  joint  pdf.  Again,  it  is  a  grid-free  algorithm  in  which  the  mass  of  fluid  is 

discretized  into  N  stochastic  particles,  the  n-th  of  these  having  position  x(n)(t)  and  composition  <J>(n)(t).  In 
each  time  step  the  composition  is  incremented  according  to  Eq.  (23),  while  the  position  is  incremented  by 

dx(n)(t)  =  2(x(n>[t]  ,t)dt  +  dxs ,  (24) 

where  the  stochastic  component  dxs  causes  a  random  walk  to  simulate  gradient  diffusion.  No 
implementations  of  this  algorithm  have  been  reported  in  the  literature. 

The  first  Monte  Carlo  algorithm  for  the  composition  joint  pdf  was  devised  by  Pope  (1981b).  In  this 
method  there  is  a  finite-difference  grid  in  physical  space.  At  each  grid  node,  the  composition  joint  pdf  is 

represented  by  N  particles,  the  n-th  having  composition  ^(n)(t).  Reaction  and  mixing  are  performed 
according  to  Eq.  (23),  while  particles  are  moved  from  node  to  node  to  simulate  convection  and  turbulent 
diffusion. 


This  algorithm  is  used  in  the  premixed  flame  calculation  of  Pope  (1981a)  and  McNutt  (1981),  and  in  the 
diffusion  flame  calculations  of  Nguyen  &  Pope  (1984),  Jones  &  Kollmann  (1987)  and  Chen  &  Kollmann 
(1988a). 


4 .  TURBULENT  FLAME  CALCULATIONS 
4.1  Turbulent  Diffusion  Flames 

Some  of  the  first  pdf  calculations  are  of  turbulent  diffusion  flames  (Frost  1975,  Janicka,  Kolbe  and 
Kollmann  1978a, b,  By  water  1982,  Nguyen  &  Pope  1984).  The  calculations  reported  by  Nguyen  &  Pope 
(1984)  arc  the  first  use  of  the  Monte  Carlo  method  for  jet  flames.  The  results  include  demonstrations  of 
convergence  of  the  solutions  as  N*1^2  tends  to  zero. 

In  the  calculations  cited  above,  the  thermochemistry  is  handled  in  a  simple  manner  —  by  assuming 
chemical  equilibrium,  for  example.  This  reduces  the  number  of  compositions  to  just  one;  namely,  the 
mixture  fraction.  Finite-rate,  multi-step  kinetics  have  been  used  by  Pope  &  Correa  (1986)  (see  also  Correa, 
Gulati  &  Pope  1988),  Jones  &  Kollmann  (1987)  and  Chen  &  Kollmann  (1988a, b).  A  computational 
challenge  is  to  implement  the  integration  of  the  rate  equation,  i.e. 


-g  =  aCsfe)  •  (25) 

in  an  efficient  manner.  All  investigators  have  used  table-look-up  algorithms. 

Considerable  attention  has  been  paid  to  the  CO/H2-air  turbulent  diffusion  flame  studied  experimentally  by 
Drake  et  al.  (1984).  Using  the  velocity-composition  joint  pdf  approach  Pope  &  Correa  (1986)  and  Correa, 
Gulati  and  Pope  (1988)  report  calculations  based  on  a  partial  equilibrium  assumption.  This  reduces  the 
number  of  compositions  to  two:  the  mixture  fraction  and  a  reaction  progress  variable  (for  the  radical 
recombination  reactions).  Again  using  the  velocity-composition  joint  pdf  approach,  Haworth,  Drake  &  Blint 
(1988)  and  Haworth,  Drake,  Pope  &  Blint  (1988)  have  made  calculations  of  this  flame  using  a  flamelet 
model. 

4.2  Turbulent  Premixed  Flames 

The  composition  joint  pdf  approach  (using  the  Monte  Carlo  method)  has  been  applied  to  premixed  flames 
by  Pope  (1981a)  and  McNutt  (1981).  The  purpose  of  the  former  calculation  was  to  demonstrate  the  ability 
of  the  pdf  method  to  handle  nonlinear  reaction  kinetics.  A  three-variable  kinetic  scheme  was  used  to  calculate 
the  oxidation  of  CO  and  formation  of  NO  in  a  propane-air  flame  stabilized  behind  a  perforated  plate. 

The  works  of  McNutt  (1981),  Pope  &  Anand  (1984)  and  Anand  &  Pope  (1987)  are  concerned  with  the 
idealized  case  of  a  statistically  steady  and  one-dimensional  turbulent  premixed  flame.  In  the  last  of  these,  the 
velocity-composition  joint  pdf  method  is  used,  and  the  effects  of  variable  density  are  studied.  It  is  shown 
that,  like  the  Bray-Moss-Libby  model  (Bray,  Libby  &  Moss  1985),  the  pdf  method  is  capable  of  accounting 
for  counter-gradient  transport  and  large  turbulence  energy  production  due  to  heat  release.  The  application  of 
the  method  to  a  spark-ignited  turbulent  flame  ball  is  described  by  Pope  &  Cheng  (1986). 

Turbulent  premixed  combustion  usually  occurs  in  the  flamelet  regime  (Pope  1987).  This  fact  presents  a 
challenge  to  any  statistical  approach,  since  the  small  scales  of  the  composition  fields  are  no  longer  governed 
by  the  turbulent  straining  motions,  but  are  determined  by  reaction  and  diffusion  occurring  in  thin  flame 
sheets.  Pope  &  Anand  (1984)  present  and  demonstrate  a  model  applicable  to  the  flamelet  regime.  But,  as 
discussed  by  Pope  (1985, 1987),  this  model  is  not  entirely  satisfactory. 

An  alternative  approach  to  treating  flamelet  combustion  is  the  stochastic  flamelet  model  of  Pope  &  Cheng 
(1988).  This  can  be  viewed  as  a  pdf  approach,  in  which  a  modelled  pdf  equation  is  solved  by  a  Monte  Carlo 


method.  In  this  case,  however,  the  pdf  is  not  that  of  fluid  properties  (i.e.  velocity  and  composition)  but  is 
rather  the  pdf  of  flamelet  properties  (i.e.  position,  area  and  orientation  of  flamelets). 

5.  DISCUSSION  AND  CONCLUSION 

The  works  reviewed  in  the  previous  section  demonstrate  that  pdf  methods  provide  a  practicable  means  of 
calculating  the  properties  of  turbulent  reactive  flows.  Calculations  have  been  made  with  thermochemical 
schemes  involving  up  to  three  composition  variables  with  finite-rate  kinetics  (e.g.  Pope  1981a,  Chen  & 
Kollmann  1988a).  And  the  Monte  Carlo  method  used  to  solve  the  pdf  equations  has  been  implemented  for  a 
variety  of  flows  including  two-dimensional  recirculating  flows  (Anand  et  al.  1989a)  and  the  three- 
dimensional  transient  flow  in  a  spark-ignition  engine  (Haworth  &  El  Tahry  1989a, b). 

The  most  advanced  method  considered  here  is  the  velocity-composition  joint  pdf  approach.  Compared  to 
moment  closures,  this  approach  has  the  advantage  that  reaction  can  be  treated  exactly  without  approximation. 
Compared  to  the  composition  joint  pdf  approach  it  has  the  advantages  that  turbulent  transport  is  treated 
exactly,  and  that  a  separate  turbulence  model  is  not  needed  to  determine  the  Reynolds  stresses. 

A  short-coming  of  the  velocity-composition  joint  pdf  approach  is  that  it  does  not  provide  a  completely 
self-contained  model,  in  that  the  turbulence  frequency  <co>  s  <e>/k  must  be  determined  by  separate  means. 

For  example,  in  some  calculations  of  simple  free  shear  layers,  it  has  been  assumed  that  <co>  is  constant 
across  the  flow,  and  scales  with  the  mean-flow  velocity  and  length  scales  (e.g.  Haworth  &  Pope  1987a, 
Pope  &  Correa  1986).  In  more  complex  flows,  another  approach  is  to  solve  the  standard  model  equation  for 

<e>  (e.g.  Haworth  &  El  Tahry  1989a)  or,  similarly,  to  solve  a  modelled  equation  for  <co>  deduced  from 
those  for  k  and  <e>  (Anand  et  al.  1989b). 

A  natural  extension  is  to  consider  f(V.\|/.C:x.t)  —  the  joint  pdf  of  velocity,  composition  and  dissipation. 
This  is  the  probability  density  function  of  the  compound  event  (H(x,t)  =  V,  <J)(x,t)  =  \j£,  e(x,t)  =  £},  where 
e(x,t)  is  the  instantaneous  mechanical  dissipation.  Following  some  preliminary  investigations  (Pope  & 

Haworth  1986,  Pope  &  Chen  1987)  a  satisfactory  model  equation  for  ffV.y.^tx.t)  has  been  developed  (Pope 
1989b).  The  incorporation  of  dissipation  within  the  pdf  allows  more  realistic  and  accurate  modelling.  But 

more  important,  the  single  equation  for  f(V.y.^:x.t)  provides  a  completely  self-contained  model  for  turbulent 
reactive  flows. 

We  now  discuss  three  areas  in  which  progress  can  be  expected  in  the  next  five  years. 

The  first  area  is  in  the  turbulent  mixing  models.  As  discussed  in  Section  2,  the  stochastic  mixing  models 
used  lead  to  the  composition  time  series  being  discontinuous  —  clearly  contrary  to  the  physics  of  the 
problem.  Nevertheless,  in  spite  of  their  lack  of  physical  appeal,  the  stochastic  models  have  many  advantages 
over  alternative  suggestions,  and,  for  inert  mixing  their  performance  is  generally  acceptable.  But  for  reactive 
flows,  especially  in  the  flamelet  regime,  their  performance  is  highly  suspect.  We  expect  that  stochastic 
models  will  be  improved  and  refined  to  account  better  for  the  microstructure  of  the  composition  fields,  and 
also  to  allow  mixing  and  reaction  to  proceed  simultaneously  at  finite  rates. 

The  second  area  of  expected  progress  is  in  the  computational  implementation  of  complex  kinetics.  When 
the  Monte  Carlo  method  is  used  to  solve  the  joint  pdf  equation  for  an  inert  flow  involving  a  compositions, 
the  computer  time  and  storage  increases  at  most  linearly  with  c.  But  with  reaction,  in  a  naive 

implementation,  for  each  particle,  on  each  time  step,  it  is  necessary  to  solve  the  coupled  set  of  a  ordinary 
differential  equations 


d0a 

“5"  =  S(x(<t>i,<|)2».*M<t>a) .  «  =  . 


(26) 


The  right-hand  side  (which  is  a  combination  of  reaction  rates)  is  computationally  expensive  to  evaluate,  and, 
as  is  well  known,  the  set  of  equations  is  likely  to  be  stiff.  Hence  such  a  naive  implementation  is 

impracticable  for  all  but  the  lowest  values  of  a. 

As  mentioned  in  Section  4.1,  the  alternative  approach  followed  by  all  investigators  is  to  implement  Eq. 
(26)  more  efficiently  through  a  different  table  look-up  scheme  (e.g.  Pope  &  Correa  1986,  Jones  &  Kollmann 
1987).  To  date  this  has  been  done  on  an  ad  hoc  basis.  It  is  expected  that  progress  will  be  made  in  the 
development  of  a  general  methodology. 

The  third  area  of  expected  progress  is  in  the  determination  of  the  mean  pressure  field  <p(x,t)>  within  the 
Monte  Carlo  algorithm.  For  thin  shear  flows,  the  mean  pressure  is  readily  determined  by  invoking  the 
boundary-layer  approximations.  For  statistically-stationary,  constant-density,  two-dimensional  recirculating 
flow,  an  algorithm  to  determine  <p>  has  been  developed  and  demonstrated  (Anand  et  al.  1989a).  But  for  the 
general  case  a  computationally  efficient  and  robust  algorithm  needs  to  be  developed.  (In  the  three- 
dimensional  transient  calculations  of  Haworth  &  El  Tahry  1989a,b,  the  Monte  Carlo  method  is  coupled  to  a 
finite-volume  code  that  determines  <p>.) 
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